Show code
library(GeoPressureR)
library(leaflet)
library(leaflet.extras)
library(raster)
library(dplyr)
library(ggplot2)
library(kableExtra)
library(plotly)
library(GeoLocTools)
setupGeolocation()
knitr::opts_chunk$set(echo = FALSE)
load(paste0("../data/1_pressure/", params$gdl_id, "_pressure_prob.Rdata"))
load(paste0("../data/2_light/", params$gdl_id, "_light_prob.Rdata"))
load(paste0("../data/3_static/", params$gdl_id, "_static_prob.Rdata"))
load(paste0("../data/4_basic_graph/", params$gdl_id, "_basic_graph.Rdata"))
Settings used
All the results produced here are generated with (1) the raw geolocator data, (2) the labeled files of pressure and light and (3) the parameters listed below.
Show code
kable(gpr) %>% scroll_box(width = "100%")
| gdl_id | crop_start | crop_end | thr_dur | extent_N | extent_W | extent_S | extent_E | map_scale | map_max_sample | map_margin | prob_map_s | prob_map_thr | shift_k | kernel_adjust | calib_lon | calib_lat | calib_1_start | calib_1_end | calib_2_start | calib_2_end | calib_2_lon | calib_2_lat | prob_light_w | thr_prob_percentile | thr_gs | RingNo | scientific_name | common_name | mass | wing_span | Color |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 16LO | 2017-01-10 | 2017-12-23 | 6 | 17 | 9 | -25 | 38 | 5 | 300 | 30 | 1.2 | 0.9 | 0 | 1.4 | 28.77132 | -22.72004 | 2017-01-10 | 2017-03-26 | 2017-11-15 | 2017-12-23 | NA | NA | 0.1 | 0.95 | 120 | NA | Halcyon senegaloides | Woodland Kingfisher | NA | NA | #FF70A6 |
Pressure timeserie
The labeling of pressure data is illustrated with this figure. The black dots indicates the pressure datapoint not considered in the matching. Each stationary period is illustrated by a different colored line.
Show code
pressure_na <- pam$pressure %>%
mutate(obs = ifelse(isoutliar | sta_id == 0, NA, obs))
p <- ggplot() +
geom_line(data = pam$pressure, aes(x = date, y = obs), colour = "grey") +
geom_point(data = subset(pam$pressure, isoutliar), aes(x = date, y = obs), colour = "black") +
# geom_line(data = pressure_na, aes(x = date, y = obs, color = factor(sta_id)), size = 0.5) +
geom_line(data = do.call("rbind", shortest_path_timeserie) %>% filter(sta_id > 0), aes(x = date, y = pressure0, col = factor(sta_id))) +
theme_bw() +
scale_colour_manual(values = pam$sta$col) +
scale_y_continuous(name = "Pressure(hPa)")
ggplotly(p, dynamicTicks = T) %>% layout(showlegend = F)
Pressure calibration
Show code
sp_pressure = do.call("rbind", shortest_path_timeserie) %>% filter(sta_id > 0)
sta_plot <- which(difftime(pam$sta$end,pam$sta$start,unit="days")>3)
par(mfrow=c(2,3))
for (i in seq_len(length(sta_plot))){
i_s = sta_plot[i]
pressure_s = pam$pressure %>%
filter(sta_id==i_s & !isoutliar)
err <- pressure_s %>% left_join(sp_pressure, by="date") %>%
mutate(
err = obs-pressure-mean(obs-pressure)
) %>% .$err
hist(err, freq = F, main = paste0("sta_id=",i_s, " | ",nrow(pressure_s)," dtpts | std=",round(sd(err),2)))
xfit <- seq(min(err), max(err), length = 40)
yfit <- dnorm(xfit, mean = mean(err), sd = sd(err))
lines(xfit, yfit, col = "red")
}
Light
Show code
raw_geolight <- pam$light %>%
transmute(
Date = date,
Light = obs
)
lightImage(tagdata = raw_geolight, offset = 0)
tsimagePoints(twl$twilight,
offset = 0, pch = 16, cex = 1.2,
col = ifelse(twl$deleted, "grey20", ifelse(twl$rise, "firebrick", "cornflowerblue"))
)
abline(v = gpr$calib_2_start, lty = 1, col = "firebrick", lwd = 1.5)
abline(v = gpr$calib_1_start, lty = 1, col = "firebrick", lwd = 1.5)
abline(v = gpr$calib_2_end, lty = 2, col = "firebrick", lwd = 1.5)
abline(v = gpr$calib_1_end, lty = 2, col = "firebrick", lwd = 1.5)
The probability map resulting from light data alone can be seen below.
Show code
li_s <- list()
l <- leaflet(width = "100%") %>%
addProviderTiles(providers$Stamen.TerrainBackground) %>%
addFullscreenControl()
for (i_r in seq_len(length(light_prob))) {
i_s <- metadata(light_prob[[i_r]])$sta_id
info <- pam$sta[pam$sta$sta_id == i_s, ]
info_str <- paste0(i_s, " | ", info$start, "->", info$end)
li_s <- append(li_s, info_str)
l <- l %>% addRasterImage(light_prob[[i_r]], opacity = 0.8, colors = "OrRd", group = info_str)
}
l %>%
addCircles(lng = gpr$calib_lon, lat = gpr$calib_lat, color = "black", opacity = 1) %>%
addLayersControl(
overlayGroups = li_s,
options = layersControlOptions(collapsed = FALSE)
) %>%
hideGroup(tail(li_s, length(li_s) - 1))
Light vs Pressure
We can compare light and pressure location at long stationary stopover (>5 days). By assuming the best match of the pressure to be the truth, we can plot the histogram of the zenith angle and compare to the fit of kernel density at the calibration site.
Show code
raw_geolight <- pam$light %>%
transmute(
Date = date,
Light = obs
)
dur <- unlist(lapply(pressure_prob, function(x) difftime(metadata(x)$temporal_extent[2],metadata(x)$temporal_extent[1], units = "days" )))
long_id <- which(dur>5)
par(mfrow = c(2, 3))
for (i_s in long_id){
twl_fl <- twl %>%
filter(!deleted) %>%
filter(twilight>shortest_path_timeserie[[i_s]]$date[1] & twilight<tail(shortest_path_timeserie[[i_s]]$date,1))
sun <- solar(twl_fl$twilight)
z_i <- refracted(zenith(sun, shortest_path_timeserie[[i_s]]$lon[1], shortest_path_timeserie[[i_s]]$lat[1]))
hist(z_i, freq = F, main = paste0("sta_id=",i_s, " | ",nrow(twl_fl),"twls"))
lines(fit_z, col = "red")
xlab("Zenith angle")
}
Similarly, we can plot the line of sunrise/sunset at the best match of pressure (yellow line) and compare to the raw and labeled light data.
Show code
lightImage(
tagdata = raw_geolight,
offset = gpr$shift_k / 60 / 60
)
tsimagePoints(twl$twilight,
offset = gpr$shift_k / 60 / 60, pch = 16, cex = 1.2,
col = ifelse(twl$deleted, "grey20", ifelse(twl$rise, "firebrick", "cornflowerblue"))
)
for (ts in shortest_path_timeserie){
if (!is.null(ts)){
twl_fl <- twl %>%
filter(twilight>ts$date[1] & twilight<tail(ts$date,1))
if (nrow(twl_fl)>0){
tsimageDeploymentLines(twl_fl$twilight,
lon = ts$lon[1], ts$lat[1],
offset = gpr$shift_k / 60 / 60, lwd = 3,col = adjustcolor("orange", alpha.f = 0.5))
}
}
}
GeoPressureViz
To visualize the path on GeoPressureViz, you will need to also load the pressure and light probability map and align them first with the code below.
Show code
sta_marginal <- unlist(lapply(static_prob_marginal, function(x) raster::metadata(x)$sta_id))
sta_pres <- unlist(lapply(pressure_prob, function(x) raster::metadata(x)$sta_id))
sta_light <- unlist(lapply(light_prob, function(x) raster::metadata(x)$sta_id))
pressure_prob <- pressure_prob[sta_pres %in% sta_marginal]
light_prob <- light_prob[sta_light %in% sta_marginal]
The code below will open with the shortest path computed with the graph approach.
Show code
geopressureviz <- list(
pam_data = pam,
static_prob = static_prob,
static_prob_marginal = static_prob_marginal,
pressure_prob = pressure_prob,
light_prob = light_prob,
pressure_timeserie = shortest_path_timeserie
)
save(geopressureviz, file = "~/geopressureviz.RData")
shiny::runApp(system.file("geopressureviz", package = "GeoPressureR"),
launch.browser = getOption("browser")
)
Stationay period information
| start | end | sta_id | col | duration |
|---|---|---|---|---|
| 2017-01-10 00:00:00 | 2017-03-25 17:20:00 | 1 | #1B9E77 | 74.7222222 days |
| 2017-03-26 03:45:00 | 2017-03-26 18:30:00 | 2 | #D95F02 | 0.6145833 days |
| 2017-03-26 20:25:00 | 2017-03-28 20:30:00 | 3 | #7570B3 | 2.0034722 days |
| 2017-03-29 03:20:00 | 2017-03-29 19:55:00 | 4 | #E7298A | 0.6909722 days |
| 2017-03-29 20:45:00 | 2017-03-30 21:45:00 | 5 | #66A61E | 1.0416667 days |
| 2017-03-31 00:25:00 | 2017-03-31 21:40:00 | 6 | #E6AB02 | 0.8854167 days |
| 2017-04-01 00:00:00 | 2017-04-01 17:50:00 | 7 | #A6761D | 0.7430556 days |
| 2017-04-01 20:35:00 | 2017-04-02 00:20:00 | 8 | #666666 | 0.1562500 days |
| 2017-04-02 01:00:00 | 2017-04-02 19:50:00 | 9 | #1B9E77 | 0.7847222 days |
| 2017-04-02 21:45:00 | 2017-04-04 01:45:00 | 10 | #D95F02 | 1.1666667 days |
| 2017-04-04 03:00:00 | 2017-04-05 01:15:00 | 11 | #7570B3 | 0.9270833 days |
| 2017-04-05 02:15:00 | 2017-04-06 23:40:00 | 12 | #E7298A | 1.8923611 days |
| 2017-04-07 01:25:00 | 2017-04-08 23:00:00 | 13 | #66A61E | 1.8993056 days |
| 2017-04-09 00:15:00 | 2017-04-09 22:15:00 | 14 | #E6AB02 | 0.9166667 days |
| 2017-04-09 22:55:00 | 2017-04-15 19:45:00 | 15 | #A6761D | 5.8680556 days |
| 2017-04-16 01:20:00 | 2017-04-21 17:55:00 | 16 | #666666 | 5.6909722 days |
| 2017-04-22 03:45:00 | 2017-04-22 21:35:00 | 17 | #1B9E77 | 0.7430556 days |
| 2017-04-23 03:45:00 | 2017-04-23 18:30:00 | 18 | #D95F02 | 0.6145833 days |
| 2017-04-24 02:05:00 | 2017-04-24 20:40:00 | 19 | #7570B3 | 0.7743056 days |
| 2017-04-24 21:20:00 | 2017-04-25 21:55:00 | 20 | #E7298A | 1.0243056 days |
| 2017-04-26 01:15:00 | 2017-04-27 00:40:00 | 21 | #66A61E | 0.9756944 days |
| 2017-04-27 03:25:00 | 2017-04-27 18:55:00 | 22 | #E6AB02 | 0.6458333 days |
| 2017-04-27 22:25:00 | 2017-04-29 00:50:00 | 23 | #A6761D | 1.1006944 days |
| 2017-04-29 02:40:00 | 2017-04-29 22:30:00 | 24 | #666666 | 0.8263889 days |
| 2017-04-30 00:30:00 | 2017-05-30 20:20:00 | 25 | #1B9E77 | 30.8263889 days |
| 2017-05-30 22:35:00 | 2017-06-01 20:25:00 | 26 | #D95F02 | 1.9097222 days |
| 2017-06-01 22:20:00 | 2017-10-16 17:10:00 | 27 | #7570B3 | 136.7847222 days |
| 2017-10-17 03:25:00 | 2017-10-17 16:30:00 | 28 | #E7298A | 0.5451389 days |
| 2017-10-18 03:25:00 | 2017-10-18 16:25:00 | 29 | #66A61E | 0.5416667 days |
| 2017-10-19 01:40:00 | 2017-10-19 18:40:00 | 30 | #E6AB02 | 0.7083333 days |
| 2017-10-19 22:40:00 | 2017-10-20 03:05:00 | 31 | #A6761D | 0.1840278 days |
| 2017-10-20 03:25:00 | 2017-10-20 18:30:00 | 32 | #666666 | 0.6284722 days |
| 2017-10-21 03:25:00 | 2017-10-21 19:15:00 | 33 | #1B9E77 | 0.6597222 days |
| 2017-10-21 23:40:00 | 2017-10-22 21:35:00 | 34 | #D95F02 | 0.9131944 days |
| 2017-10-23 00:30:00 | 2017-10-24 00:05:00 | 35 | #7570B3 | 0.9826389 days |
| 2017-10-24 02:25:00 | 2017-10-25 00:20:00 | 36 | #E7298A | 0.9131944 days |
| 2017-10-25 02:25:00 | 2017-10-26 18:40:00 | 37 | #66A61E | 1.6770833 days |
| 2017-10-26 21:00:00 | 2017-10-27 21:25:00 | 38 | #E6AB02 | 1.0173611 days |
| 2017-10-28 01:20:00 | 2017-10-31 18:45:00 | 39 | #A6761D | 3.7256944 days |
| 2017-10-31 23:30:00 | 2017-11-01 17:30:00 | 40 | #666666 | 0.7500000 days |
| 2017-11-01 23:55:00 | 2017-11-02 19:45:00 | 41 | #1B9E77 | 0.8263889 days |
| 2017-11-02 22:25:00 | 2017-11-02 23:55:00 | 42 | #D95F02 | 0.0625000 days |
| 2017-11-03 02:10:00 | 2017-11-04 01:40:00 | 43 | #7570B3 | 0.9791667 days |
| 2017-11-04 02:50:00 | 2017-11-07 23:10:00 | 44 | #E7298A | 3.8472222 days |
| 2017-11-08 01:00:00 | 2017-11-08 23:10:00 | 45 | #66A61E | 0.9236111 days |
| 2017-11-09 00:15:00 | 2017-11-09 19:40:00 | 46 | #E6AB02 | 0.8090278 days |
| 2017-11-09 21:00:00 | 2017-11-10 21:25:00 | 47 | #A6761D | 1.0173611 days |
| 2017-11-11 00:10:00 | 2017-11-12 17:50:00 | 48 | #666666 | 1.7361111 days |
| 2017-11-12 22:15:00 | 2017-11-13 20:20:00 | 49 | #1B9E77 | 0.9201389 days |
| 2017-11-14 01:30:00 | 2017-11-14 19:00:00 | 50 | #D95F02 | 0.7291667 days |
| 2017-11-15 00:15:00 | 2017-12-22 23:55:00 | 51 | #7570B3 | 37.9861111 days |